{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Feature Engineering Notebook Three: Dimensionality Reduction**  \n",
    "*Author: Yingxiang Chen, Zihan Yang*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Reference**\n",
    "- https://scikit-learn.org/stable/modules/decomposition.html#principal-component-analysis-pca\n",
    "- https://sebastianraschka.com/faq/docs/lda-vs-pca.html\n",
    "- https://en.wikipedia.org/wiki/Linear_discriminant_analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Dimension-Reduction\" data-toc-modified-id=\"Dimension-Reduction-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Dimension Reduction</a></span><ul class=\"toc-item\"><li><span><a href=\"#Unsupervised-Methods\" data-toc-modified-id=\"Unsupervised-Methods-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Unsupervised Methods</a></span><ul class=\"toc-item\"><li><span><a href=\"#PCA-(Principal-Components-Analysis)\" data-toc-modified-id=\"PCA-(Principal-Components-Analysis)-1.1.1\"><span class=\"toc-item-num\">1.1.1&nbsp;&nbsp;</span>PCA (Principal Components Analysis)</a></span></li></ul></li><li><span><a href=\"#Supervised-Methods\" data-toc-modified-id=\"Supervised-Methods-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Supervised Methods</a></span><ul class=\"toc-item\"><li><span><a href=\"#LDA-(Linear-Discriminant-Analysis)\" data-toc-modified-id=\"LDA-(Linear-Discriminant-Analysis)-1.2.1\"><span class=\"toc-item-num\">1.2.1&nbsp;&nbsp;</span>LDA (Linear Discriminant Analysis)</a></span></li></ul></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimension Reduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After data preprocessing and feature selection, we have generated a good feature subset. But sometimes, this subset might still contain too many features and cost so much computing power to train. In this case, we can use dimension reduction techniques to further compress our feature subset. But this might deprecate model performance.\n",
    "\n",
    "We can also apply dimension reduction methods directly after data preprocessing if we don't have much time on feature selection. The dimension reduction algorithm can compress the original feature space and generate a feature subset for us.\n",
    "\n",
    "Specifically, we will introduce PCA and LDA (Linear Discriminant Analysis)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unsupervised Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCA (Principal Components Analysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA is an **unsupervised** technique that finds the directions of maximal variance. It uses a few unrelated features to represent original features in the dataset and tries to retain as much information (variance) as possible. More math detail can be viewed from a [repo](https://github.com/YC-Coder-Chen/Unsupervised-Notes/blob/master/PCA.md) written by us in Github."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-30T03:37:18.226607Z",
     "start_time": "2020-03-30T03:37:15.126399Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "# load dataset\n",
    "from sklearn.datasets import fetch_california_housing\n",
    "dataset = fetch_california_housing()\n",
    "X, y = dataset.data, dataset.target # use california_housing dataset as example\n",
    "\n",
    "# use the first 20000 observations as train_set\n",
    "# the rest observations as test_set\n",
    "train_set = X[0:20000,]\n",
    "test_set = X[20000:,]\n",
    "train_y = y[0:20000]\n",
    "\n",
    "# we need to standardize the data first or the PCA will only comopress features in\n",
    "# large scale\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "model = StandardScaler()\n",
    "model.fit(train_set) \n",
    "standardized_train = model.transform(train_set)\n",
    "standardized_test = model.transform(test_set)\n",
    "\n",
    "# start compressing\n",
    "compressor = PCA(n_components=0.9) # set n_components=0.9 =>\n",
    "# select the number of components such that the amount of variance\n",
    "# explained is greater than 90% of the original variance\n",
    "# we can also set n_components to be the number of features we want directly\n",
    "\n",
    "compressor.fit(standardized_train) # fit on trainset\n",
    "transformed_trainset = compressor.transform(standardized_train) # transform trainset (20000,5)\n",
    "transformed_testset = compressor.transform(standardized_test) # transform test set\n",
    "\n",
    "assert transformed_trainset.shape[1] == transformed_testset.shape[1] # same number of features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-30T03:37:18.780647Z",
     "start_time": "2020-03-30T03:37:18.242698Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8FeXZ//HPxSo7su8EIYDsYAQV6opWrXWviisuRdtStba/51Hroy192sdareIuUtwrdS9aLSiKigsQVFYFAwSJoCD7FrJdvz9mkh5jSCbI5CQ53/frlRdnzrln5joTcq4zM/d93ebuiIiIANRJdgAiIlJ9KCmIiEgJJQURESmhpCAiIiWUFEREpISSgoiIlFBSEBGREkoKIiJSQklBRERK1Et2AJXVpk0bT0tLS3YYIiI1yvz5879x97YVtatxSSEtLY3MzMxkhyEiUqOY2eoo7XT5SERESigpiIhICSUFEREpoaQgIiIllBRERKRErEnBzE40s2VmlmVm15fxenczm2lmC81slpl1iTMeEREpX2xJwczqAvcBJwH9gDFm1q9Us9uBx919EDAB+L+44hERkYrFeaYwHMhy95XungdMBU4r1aYfMDN8/FYZr4uIpLzNO/P4y/TPyP5mZ+z7inPwWmdgTcJyDjCiVJsFwFnAROAMoJmZtXb3jTHGJSJSI2zZlcfkd1fx6PvZ7MwroEOLRqS1aRLrPuNMClbGc15q+TfAvWY2FngH+BIo+M6GzMYB4wC6deu2f6MUEalmtu7KZ/LslTzyXjY79hTwo4Edufq4dPp0aBb7vuNMCjlA14TlLsDaxAbuvhY4E8DMmgJnufvW0hty90nAJICMjIzSiUVEpFbYuiufv723ikdmr2L7ngJOHtiBq49Lp2+H5lUWQ5xJYR6QbmY9CM4AzgPOT2xgZm2ATe5eBNwATIkxHhGRamnr7nymzF7FlPdWsT23gJMGBMng4I5VlwyKxZYU3L3AzMYD04G6wBR3X2JmE4BMd58GHA38n5k5weWjX8QVj4hIdbMtN0gGf5sdJIMf9m/PNcf1pl+nqk8Gxcy9Zl2NycjIcFVJFZGabFtuPo++l83kd1eyLbeAE/q155rR6fTv1CK2fZrZfHfPqKhdjSudLSJSU20vTgazV7F1dz7H92vPNcelM6BzfMmgspQURERitj03n8fez+bhd4NkMPrgdlw7une1SgbFlBRERGKyY09BmAxWsmVXPsf1DZLBwC7VLxkUU1IQEdnPdu4p4LEPsnn4nZVs3pXPsX3bcc1x6Qzu2jLZoVVISUFEZD/ZuaeAxz9YzaR3VrB5Vz5H92nLtaN7M6QGJINiSgoiIt/TrrziZLCSTTvzOKp3W64dnc7QbgcmO7RKU1IQEdlHu/IKePLD1Tz09ko27szjyN5tuea4dA7pXvOSQTElBRGRStqdVxgkg3dW8M2OPH6Q3oZrR6dzSPdWyQ7te1NSEBGJaHdeIU/NWc2DbwfJYFSvIBlkpNX8ZFBMSUFEpAK5+YU8NecLHnx7BRu272Fkr9Y8MLo3h9aiZFBMSUFEZC9y8wv5+5wveCBMBkf0bM195w9jeI/alwyKKSmIiJSSm1/I03O/4IFZK1i/fQ+HHdSKe8YM5bCDWic7tNgpKYiIhHLzC/nHvDXcPyuLr7ftYXiPVkw8byiH96z9yaCYkoKIpLw9BWEyeGsFX23LZXhaK+48dwhH9GyT7NCqnJKCiKSsPQWFPDNvDffPWsG6rbkcmnYgfz1nMIf3bI1ZWTMK135KCiKScvYUFPJsZg73v5XF2q25ZHQ/kNt/MpgjUjgZFFNSEJGUkVdQxLPz13Dfm0EyGNatJX8+exCjerVJ+WRQLFJSMLNRQLq7P2JmbYGm7r4q3tBERPaPvIIinpufw31vZfHllt0M7daSW88axA/SlQxKqzApmNktQAbQB3gEqA88CYyMNzQRke8nv7CI5+fncM+bQTIY0rUlfzpzIEcqGexVlDOFM4ChwEcA7r7WzJrFGpWIyPeQX1jECx8FySBn824Gd23J/54xgKN7t1UyqECUpJDn7m5mDmBmTWKOSURkn+QXFvHiR19yz1ufs2bTbgZ1acEfThvA0X2UDKKKkhSeMbOHgJZm9lPgMuDheMMSEYmuoLCIFz/+knvezOKLTbsY2LkFvx/bn2P6tFMyqKQKk4K7325mxwPbCO4r3Ozur0fZuJmdCEwE6gKT3f3WUq93Ax4DWoZtrnf3Vyv3FkQkVRUUFvHSJ2u5583PWb1xFwM6N2fyxRkcd7CSwb6KcqO5B/BucSIws0Zmlubu2RWsVxe4DzgeyAHmmdk0d1+a0Owm4Bl3f8DM+gGvAmn79E5EJGUUFBbxzzAZZG/cRf9OzXn44gxGKxl8b1EuHz0LHJGwXBg+d2gF6w0Hstx9JYCZTQVOAxKTggPNw8ctgLUR4hGRFFVY5Exb8CV3z8xi1Tc76dexOZMuOoTj+7VXMthPoiSFeu6eV7zg7nlm1iDCep2BNQnLOcCIUm1+B8wws18CTYDRZW3IzMYB4wC6desWYdciUpsUFjkvL1jL3TM/Z+U3Ozm4Y3MeuugQTlAy2O+iJIUNZnaqu08DMLPTgG8irFfWb8pLLY8BHnX3O8zscOAJMxvg7kXfWsl9EjAJICMjo/Q2RKSWKixyXlm4lokzP2flhp307dCMBy8cxgn9OlCnjpJBHKIkhauAp8zsXoIP+jXAxRHWywG6Jix34buXhy4HTgRw9w/M7ACgDbA+wvZFpJYqLHL+tWgdd8/8nKz1O+jTvhkPXDCMH/ZXMohblN5HK4DDzKwpYO6+PeK25wHp4Y3qL4HzgPNLtfkCOA541MwOBg4ANkQNXkRql6KEZPD5+h30bt+U+y8YxolKBlUmSu+jhsBZBL2C6hVfv3P3CeWt5+4FZjYemE7Q3XSKuy8xswlAZng56tfAw2b2K4JLS2PdXZeHRFJMUZHz6uJ1THwjSAbp7Zpy7/lDOXlARyWDKhbl8tE/ga3AfGBPZTYejjl4tdRzNyc8XopqKImkrKIi57XFXzFx5nKWf72DXu2acs+YofxooJJBskRJCl3c/cTYIxGRlFFU5Exf8hUTZ37OZ19tp2fbJtwdJoO6SgZJFSUpvG9mA919UezRiEit5u68tWw9d8xYzpK12ziobRMmnjeEUwZ1UjKoJqIkhVHAWDNbRXD5yAB390GxRiYitcp7Wd9w+4xlfPzFFrq1asxfzxnMaUM6KxlUM1GSwkmxRyEitdb81Zu4ffpyPli5kY4tDuD/zhzI2Yd0oX7dOskOTcoQpUvqagAza0fQZVREpEKLv9zK7TOWMWvZBto0bcgtP+7HmOHdOKB+3WSHJuWI0iX1VOAOoBPBoLLuwKdA/3hDE5GaaPnX2/nrjOX8e8lXtGhUn/8+sS+XHNGdxg00JXxNEOW39AfgMOANdx9qZscQlKcQESmR/c1O7npjOf9csJYmDepx7eh0LhvVg+YH1E92aFIJUZJCvrtvNLM6ZlbH3d8ysz/HHpmI1AhfbtnNPTM/59n5OdSva1x5ZE+uPPIgDmwSpW6mVDdRksKWsMTFOwQ1kNYDBfGGJSLV3fptudz3VhZPzw2KIV90WHd+fkxP2jXTrceaLEpSOA3IBX4FXEAw70G5JS5EpPbavDOPB99ewWMfZFNQ6Pwkoyu/PLYXnVo2SnZosh9E6X20M2HxsRhjEZFqbFtuPpPfXcWU2avYmVfAGUM6c83odLq3bpLs0GQ/2mtSMLPZ7j7KzLbz7XkQigevNd/LqiJSi+zKK+DR97N56O2VbN2dz8kDO/Cr0b1Jb98s2aFJDPaaFNx9VPivfvMiKSg3v5Cn5nzBA7Oy+GZHHsf2bcd1x/dmQOcWyQ5NYlTu5SMzqwMsdPcBVRSPiCRZXkERz85fwz0zs/hqWy4je7XmoeP7cEj3A5MdmlSBcpOCuxeZ2QIz6+buX1RVUCJS9QqLnJc+/pK7Zi5nzabdDOvWkr+eO5gjerZJdmhShaL0PuoILDGzuUDJTWd3PzW2qESkyhRPcHPn68tZsWEn/Ts155GxAzi6T1uKJ9WS1BElKfw+9ihEpMq5OzM/Xc8dry/n03XbSG/XlAcvDOZBVjJIXVG6pL5dFYGISNVwd97L2sjtM5bxyZotdG/dmLvOHcKPB2tOA4lWEO8w4B7gYKABwXzLO9UlVaTmyczexF+mL2POqk10anEAt545kLNUxloSRLl8dC9wHvAskAFcDKTHGZSI7F+LcoIy1m8v30DbZg35/an9OW94VxrWUxlr+bZItWzdPcvM6rp7IfCImb0fc1wish8s+2o7f319GdOXfE3LxvW54aS+XHx4Go0aKBlI2aIkhV1m1gD4xMxuA9YBkca1m9mJwESCS06T3f3WUq/fCRwTLjYG2rl7y6jBi0jZVoVlrKctWEvTBvX41ejeXDYqjWYqYy0ViJIULgLqAOMJiuJ1Bc6qaCUzqwvcBxwP5ADzzGyauy8tbuPuv0po/0tgaKWiF5Fvydm8i3tmZvHcRzk0qFuHq44Kyli3bKwy1hJNlKQwDHjV3bdRue6pw4Esd18JYGZTCSquLt1L+zHALZXYvoiEvi4pY/0FZsYlh6fxs6N70rZZw2SHJjVMlKRwKnCXmb0DTAWmu3uU+RQ6A2sSlnOAEWU1NLPuQA/gzQjbFZHQpuIy1u9nU1jknHNoV8YfozLWsu+ijFO41MzqAycB5wP3m9nr7n5FBauW1eHZy3gOgt5Nz4U3sr+7IbNxwDiAbt26VRSySK23dXc+k99dyZTZq9idX8jpQztzzXEqYy3fX9TeR/lm9hrBh3ojgstAFSWFHIL7D8W6AGv30vY84Bfl7H8SMAkgIyNjb4lFpNbbuae4jPUKtuUW8KOBHfnV8en0aqdixrJ/RBm8diLBh/YxwCxgMnBOhG3PA9LNrAfwZbiN88vYfh/gQOCDyFGLpJjc/EKe/HA1D8xawcadeRzXtx3XndCb/p1Uxlr2ryhnCmMJ7iVc6e57om7Y3QvMbDwwnaBL6hR3X2JmE4BMd58WNh0DTHV3nQGIlJJXUMQzmWu4982gjPWoXm247oTeDOumMtYSD6tpn8UZGRmemZmZ7DBEYlVQWMRLn6xlYljGOqP7gfz6hD4c3rN1skOTGsrM5rt7RkXtIt1TEJGqUVTk/GvROu58YzkrN+xkYOcW/OHSARzVW2WspWooKYhUA+7OG5+u544Zy/jsq+30bt+UBy88hB/2b69kIFVKSUEkidyd2VnfcPuM5SxYs4W01o2ZeN4QThmkMtaSHHtNCma2iL2PK8DdB8USkUiKmLtqE7fPWMbcVZvo3LIRt501iDOHdaaeylhLEpV3pnBK+G/x+IEnwn8vAHbFFpFILbdgzRbueH0574RlrCec1p9zD1UZa6ke9poU3H01gJmNdPeRCS9db2bvARPiDk6kNvnsq23cMWM5ry/9mgMb1+fGk/ty0WEqYy3VS5R7Ck3MbJS7zwYwsyOIWDpbRGDlhh3c+cbnvLIwKGN93fG9uWxUD5o21C09qX6i/K+8HJhiZi0I7jFsBS6LNSqRWmDNpl3cPfNznv8oh4b16vKzo3oyTmWspZqLUhBvPjDYzJoTDHbbGn9YIjXX19tyuffNLKbOC8pYXzqyBz87uidtmqqMtVR/UWoftQf+BHRy95PMrB9wuLv/LfboRGqQjTv28ODbK3j8g9UUFjnnHtqV8cf2omMLlbGWmiPK5aNHgUeA34bLy4F/AEoKIgQlKR59P5s7X1/O7vxCzhjahWtHp9O1VeNkhyZSaVGSQht3f8bMboCSQndlznsgkmqWrN3K9c8vYtGXWzm2bztuPPlgerVrmuywRPZZlKSw08xaEw5kM7PDCG42i6Ss3XmF3DVzOZPfXcWBjetzz5ihnDKoo0pSSI0XJSlcB0wDeobjE9oCZ8calUg1Nvvzb7jxxUV8sWkX52R04caTD1aPIqk1ovQ++sjMjgL6EEyxuczd82OPTKSa2bwzj//916c8/1EOaa0b8/efjuCInm2SHZbIfhV19MxwIC1sP8zMcPfHY4tKpBpxd6YtWMuEl5eydXc+Pz+6J1cfl84B9TUSWWqfKF1SnwB6Ap8AxTeYHVBSkFovZ/MubnppMbOWbWBwlxY8ecUIDu7YPNlhicQmyplCBtBP02VKKiksch59P5s7ZiwD4OZT+nHJEWkqZy21XpSksBjoAKyLORaRamHp2m3c8MJCFuRs5Zg+bfnD6QPocqDGHEhqiDROAVhqZnOBPcVPuvupsUUlkgS5+YVMnPk5k95ZSctG9bl7zFB+rG6mkmKiJIXfxR2ESLK9nxV0M83euIuzD+nCb08+mAObqJuppJ4oXVLfropARJJhy648/vivT3l2fg7dWzfmqStGMLKXuplK6ipvOs7Z7j7KzLbz7Wk5DXB3r7ALhpmdCEwE6gKT3f3WMtqcQ3A24sACdz+/cm9BpPLcnZcXrmPCy0vYvCufq47qybWj1c1UpLyZ10aF/zbblw2bWV3gPuB4IAeYZ2bT3H1pQpt04AZgpLtvNrN2+7Ivkcr4cstubnpxEW8t28CgLi147LLh9O/UItlhiVQLkad+Cj+wDyhedvcvKlhlOJDl7ivD9acCpwFLE9r8FLjP3TeH21wfNR6Ryiosch57P5vbZyzDHW760cFcOrKHupmKJIgyeO1U4A6gE7Ae6A58CvSvYNXOwJqE5RxgRKk2vcN9vEdwiel37v7vMmIYB4wD6NatW0Uhi3zHZ19t47+fX8SCNVs4qndb/vf0ASptLVKGKGcKfwAOA95w96FmdgwwJsJ6ZX39Kj0Arh6QDhwNdAHeNbMB7r7lWyu5TwImAWRkZGgQnUSWm1/IPW9+zkNvr6RFo/pMPG8Ipw7upG6mInsRJSnku/tGM6tjZnXc/S0z+3OE9XKArgnLXYC1ZbT5MCywt8rMlhEkiXlRghcpzwcrNnLji4tY9c1OzhrWhZt+pG6mIhWJkhS2mFlT4B3gKTNbDxREWG8ekG5mPYAvgfOA0j2LXiI463jUzNoQXE5aGTV4kbJs3ZXPn179lH9krqFbq8Y8efkIRqWrm6lIFFGSwmlALvAr4AKgBTChopXCGdrGA9MJ7hdMcfclZjYByHT3aeFrJ5jZUoJie//P3Tfu21uRVOfu/GvROn43bSmbd+Vx5VEHce1xvWnUQN1MRaKymlbnLiMjwzMzM5MdhlQza7fs5n9eWszMz9YzoHNzbj1zEAM6q5upSDEzm+/uGRW1K2/wWpmD1qjE4DWRuBUWOU98kM1fpi+jyOG3Jx/MpSPTqFe3TrJDE6mRyhu8tk+D1kSqyrKvtvPfzy/kkzVb+EF6G/50xkB1MxX5niINXjOzYcAogjOF2e7+caxRiZQjN7+Qe9/M4sG3V9C8UX3uPHcwpw/prG6mIvtBlMFrNwM/AV4In3rUzJ519/+NNTKRMny4ciM3vrCIld/s5MyhnbnplH60UjdTkf0mypnCGGCou+cCmNmtwEeAkoJUma2787n1tU95eu4aurZqxOOXDefI3m2THZZIrRMlKWQT1DzKDZcbAiviCkgkkbvz2uKvuGXaEjbu2MO4Iw/i2tHpNG4QuWyXiFRClL+sPcASM3ud4J7C8cBsM7sbwN2vjjE+SWHrtu7mf15awhuffk3/Ts15ZOyh6mYqErMoSeHF8KfYrHhCEQkUFTlPzlnNbf9eRkFRETee3JfLRvZQN1ORKhAlKbxWuqS1mfVx92UxxSQpbPnX27n++YV89EXQzfSPpw+kW2t1MxWpKlGSwrtm9j/u/gyAmf0auBzoF2tkklJy8wu5/60sHnh7BU0b1uOOnwzmzGHqZipS1aIkhaOBSWb2E6A9wVwKw+MMSlLL3FWbuP6FhazcsJPTh3Tif07pR+umDZMdlkhKqjApuPs6M/s3wbSZRcAN7r4j9sik1gu6mX7G03O/oMuBjXjssuEcpW6mIkkVZfDa68A6YADBnAhTzOwdd/9N3MFJ7fXvxeu4+Z9L+GbHHq4Y1YPrTuitbqYi1UCUv8L73P2l8PEWMzuC4KxBpNK+2prLzf9czIylX9OvY3MmX5LBoC4tkx2WiISiXD56ycy6A+nu/gZQH7gr9sikVikqcp6a+wW3vfYZeYVFXH9SXy4f1YP66mYqUq1EuXz0U2Ac0AroSXAJ6UHguHhDk9ri86+3c8MLi8hcvZmRvVrzpzMG0r11k2SHJSJliHL56BcEvY3mALj752bWLtaopFbYU1DI/W+t4P5ZWTRpWI/bfzKYs9TNVKRai1Tmwt3ziv+Qzawe3558R+Q75mVv4oYXFpG1fgenhd1M26ibqUi1FyUpvG1mNwKNzOx44OfAy/GGJTXVttx8/vzaZzw15ws6t2zEI5ceyjF9dGIpUlNESQrXE4xgXgRcCbwKTI4zKKmZ/r34K26ZtpgN2/dw2cge/PqE3jRpqG6mIjVJlN5HRcDD4Y/Id3y9LehmOn3J1/Tt0IxJF2UwuKu6mYrURLF+jTOzE4GJQF1gsrvfWur1scBfgC/Dp+51d52F1BBFRc7f537Bn8Nupv91Yh9++oOD1M1UpAaLLSmYWV3gPoL5F3KAeWY2zd2Xlmr6D3cfH1ccEo+s9Tu44YWFzMvezBE9g26maW3UzVSkpoucFMysibvvrMS2hwNZ7r4yXH8qcBpQOilIDZJXUMQDs1Zw31tZNGpQl9vOHsRPDumibqYitUSF5/lmdoSZLSWojoqZDTaz+yNsuzOwJmE5J3yutLPMbKGZPWdmXaMELckxf/UmfnT3u9z5xnJ+OKADb1x3FOdkdFVCEKlFopwp3An8EJgG4O4LzOzICOuV9UlRenzDy8DT7r7HzK4CHgOO/c6GzMYRjKqmW7duEXYt+9P23Hxu+/cynpyzmk4tGvHI2EM5pq+6mYrURpEuH7n7mlLfBgsjrJYDJH7z7wKsLbXdjQmLDwN/3sv+JwGTADIyMjRwrgrNWPIVN/9zCV9vz2XsEWn85oQ+6mYqUotF+eteE1ZGdTNrAFxNeCmpAvOAdDPrQdC76Dzg/MQGZtbR3deFi6dG3K5Ugc0787jxxUW8tvgr+nZoxoMXHcIQdTMVqfWiJIWrCLqVdib49j+DoB5Sudy9wMzGA9MJuqROcfclZjYByHT3acDVZnYqUABsAsbu07uQ/errbblcOHkOqzfu4v/9sA/jjlQ3U5FUYe7lX40xs7buvqGK4qlQRkaGZ2ZmJjuMWmvNpl1cMHkOG3fs4eFLMjiiZ5tkhyQi+4GZzXf3jIraRTlTeN/MVgH/AJ539y3fOzqplj7/ejsX/m0OuflFPPXTw3S5SCQFVXhNwN3TgZuA/sBHZvaKmV0Ye2RSpRblbOWchz6gyOGZKw9XQhBJUZEuFLv7XHe/jmBA2iaCrqNSS8xdtYnzH/6Qxg3q8eyVh9OnQ7NkhyQiSRJl8FpzM7vEzF4D3gfWESQHqQVmLVvPxVPm0K55Q5772eEqVSGS4qLcU1gAvARMcPcPYo5HqtCri9ZxzdSP6d2+GY9fNpzWmgRHJOVFSQoHeUVdlKTGeSZzDdc/v5Bh3Q7kb2MPpUWj+skOSUSqgb0mBTO7y92vBaaZ2XeSgrufGmtkEpsps1cx4ZWl/CC9DQ9ddAiNG2iEsogEyvs0eCL89/aqCETi5+7cPTOLO99Yzon9OzBxzBAa1qub7LBEpBrZa1Jw9/nhwyHuPjHxNTO7Bng7zsBk/3J3/vivT5k8exVnDevCn88aSD2NUhaRUqJ8KlxSxnNj93McEqPCIueGFxYxefYqxh6Rxl/OHqSEICJlKu+ewhiCAnY9zGxawkvNgI1lryXVTV5BEdc98wmvLFzHL4/txXXH99b8ByKyV+XdUygek9AGuCPh+e3AwjiDkv0jN7+Qnz05n7eWbeDGk/sy7sieyQ5JRKq58u4prAZWA4dXXTiyv2zPzeeKxzKZm72JP50xkPNHaHIiEalYlBHNh5nZPDPbYWZ5ZlZoZtuqIjjZN5t35nHB5DnMX72Zu84dooQgIpFF6aB+L8EEOc8CGcDFQK84g5J9VzIXwqZdPHTRIRx3cPtkhyQiNUjU6TizzKyuuxcCj5jZ+zHHJfsgcS6ERy89VHMhiEilRUkKu8JpOD8xs9sIbj6ralo1o7kQRGR/iNJZ/SKC6TTHAzuBrsBZcQYllaO5EERkf6nwTCHshQSwG/h9vOFIZc1dtYnLH51H80b1eeqKESp9LSLfS3mD1xYBe62O6u6DYolIIpu1bD1XPTmfzi0b8eQVI+jYolGyQxKRGq68M4VTqiwKqTTNhSAicaho8JpUQ5oLQUTiEmXw2nYz2xb+5FZm8JqZnWhmy8wsy8yuL6fd2WbmZpZRmeBT0ZTZq/iv5xYyslcbHr98uBKCiOxXUW40f2sWdzM7nQhzNJtZXeA+4HggB5hnZtPcfWmpds2Aq4E5lYg75WguBBGpCpWun+zuLwHHRmg6HMhy95XungdMBU4ro90fgNuA3MrGkiqK50K4843lnDWsC/eeP1QJQURiUeGZgpmdmbBYh6DURZQ5mzsDaxKWc4ARpbY9FOjq7q+Y2W/KiWEcMA6gW7fUquNTWOT89sVFTJ23hrFHpHHzKf2oU0elr0UkHlFGNP844XEBkE3Z3/hLK+uTqySZmFkd4E4iTNjj7pOASQAZGRlRElKtoLkQRKSqRbmncOk+bjuHYPRzsS7A2oTlZsAAYFb4QdcBmGZmp7p75j7us9bQXAgikgxRLh/1AH4JpCW2d/dTK1h1HpAerv8lQaXV8xPW30owgU/xfmYBv1FC0FwIIpI8US4fvQT8DXgZKIq6YXcvMLPxwHSC2klT3H2JmU0AMt19WvlbSE2bd+ZxySNzWbp2G3edO4TThnROdkgikkKiJIVcd797Xzbu7q8Cr5Z67ua9tD16X/ZRm2guBBFJtihJYaKZ3QLMAPYUP+nuH8UWVQrSXAgiUh1ESQoDCcpnH8t/Lh850cYqSASaC0FEqosoSeEM4KBwAJrsZ4tytnLxlDnUq1uHZ648nD4dmlW8kohITKKMaF4A6KtrDOau2sT5D39I4wa9Zf3HAAANI0lEQVT1eFYJQUSqgShnCu2Bz8xsHt++p1BRl1Qph+ZCEJHqKEpSuCX2KFKM5kIQkeoqyojmt6sikFShuRBEpDqLMqJ5O/+pWdQAqA/sdPfmcQZWG02ZvYoJryzlB+lteOiiQ2jcIMqJmohI1YltPgX5j8S5EH7Yvz13j1HpaxGpnuKcT0H49lwIZw7rzH3nD1NCEJFqK875FFJe4lwIlxzenVt+3F9zIYhItRbnfAopLXEuhPHH9OLXJ2guBBGp/uKcTyFlJc6FcMNJfbnyKM2FICI1Q4X3FMzsMTNrmbB8oJlNiTesmmt7bj6XTJnLrOUb+OMZA5QQRKRGiXL5aJC7bylecPfN4dzKUkrxXAhLNBeCiNRQUZJCHTM70N03A5hZq4jrpZRvzYVw4SGM7qe5EESk5ony4X4H8L6ZPUfQ6+gc4I+xRlXDaC4EEaktotxoftzMMgnGJhhwprsvjT2yGiJxLoQnrxjB0G4HJjskEZF9FukyUJgElAhKKZ4LoW6dOvzjysPo20GVP0SkZtO9gX00d9UmLn90Hs0b1eepK0aQ1qZJskMSEfnelBT2QfFcCJ1aNuLJy0fQqaXmQhCR2qHStY8qw8xONLNlZpZlZteX8fpVZrbIzD4xs9lm1i/OePaHVxet46ePZ3JQm6Y8c+XhSggiUqvElhTMrC5wH3AS0A8YU8aH/t/dfaC7DwFuA/4aVzz7wzOZaxj/948Y1KUlT487jDaaHEdEapk4zxSGA1nuvtLd84CplKqZ5O7bEhabUI0L7U2ZvYr/em4hI3u14YnLh2tyHBGpleK8p9AZWJOwnAOMKN3IzH4BXEcwgU+1K8mtuRBEJJXEeaZQVknQ75wJuPt97t4T+G/gpjI3ZDbOzDLNLHPDhg37Ocy901wIIpJq4kwKOUDXhOUuwNpy2k8FTi/rBXef5O4Z7p7Rtm3b/Rji3hUWOTe8sIjJs1dxyeHduf3swdSrG+t9eRGRpIvzU24ekG5mPcysAXAeMC2xgZmlJyz+CPg8xngiyyso4pqpHzN13hrGH9OL352qyXFEJDXEdk/B3QvMbDwwHagLTHH3JWY2Ach092nAeDMbDeQDm4FL4oonKs2FICKpLNbBa+7+KvBqqeduTnh8TZz7r6ztuflc8Vgmc7M38cczBnDBiO7JDklEpEppRHNIcyGIiCgpALB+Wy4X/m0O2Rs1F4KIpLaUTwqaC0FE5D9SOilkrd/OhZPnsju/UHMhiIiQwklh8ZdbuXjKXOqYaS4EEZFQSiaFedmbuOwRzYUgIlJayiWFt5dv4MonMjUXgohIGVIqKby2aB1XT/2Y9HbNePzy4Sp9LSJSSsokhRc+yuE3zy5gaLcDmTL2UJW+FhEpQ8okhW6tGjP64Pbcdd4QGjdImbctIlIpKfPpmJHWioy0VskOQ0SkWlMtaBERKaGkICIiJZQURESkhJKCiIiUUFIQEZESSgoiIlJCSUFEREooKYiISAlz92THUClmtgFYvY+rtwG+2Y/h7C+Kq3IUV+VV19gUV+V8n7i6u3vbihrVuKTwfZhZprtnJDuO0hRX5SiuyquusSmuyqmKuHT5SERESigpiIhIiVRLCpOSHcBeKK7KUVyVV11jU1yVE3tcKXVPQUREypdqZwoiIlKOWpcUzGyKma03s8V7ed3M7G4zyzKzhWY2rJrEdbSZbTWzT8Kfm6sorq5m9paZfWpmS8zsmjLaVPkxixhXlR8zMzvAzOaa2YIwrt+X0aahmf0jPF5zzCytmsQ11sw2JByvK+KOK2Hfdc3sYzN7pYzXqvx4RYwrmccr28wWhfvNLOP1+P4m3b1W/QBHAsOAxXt5/WTgNcCAw4A51SSuo4FXknC8OgLDwsfNgOVAv2Qfs4hxVfkxC49B0/BxfWAOcFipNj8HHgwfnwf8o5rENRa4t6r/j4X7vg74e1m/r2Qcr4hxJfN4ZQNtynk9tr/JWnem4O7vAJvKaXIa8LgHPgRamlnHahBXUrj7Onf/KHy8HfgU6FyqWZUfs4hxVbnwGOwIF+uHP6VvzJ0GPBY+fg44zsysGsSVFGbWBfgRMHkvTar8eEWMqzqL7W+y1iWFCDoDaxKWc6gGHzahw8PT/9fMrH9V7zw8bR9K8C0zUVKPWTlxQRKOWXjJ4RNgPfC6u+/1eLl7AbAVaF0N4gI4K7zc8JyZdY07ptBdwH8BRXt5PSnHK0JckJzjBUFCn2Fm881sXBmvx/Y3mYpJoaxvINXhG9VHBMPQBwP3AC9V5c7NrCnwPHCtu28r/XIZq1TJMasgrqQcM3cvdPchQBdguJkNKNUkKccrQlwvA2nuPgh4g/98O4+NmZ0CrHf3+eU1K+O5WI9XxLiq/HglGOnuw4CTgF+Y2ZGlXo/tmKViUsgBEjN+F2BtkmIp4e7bik//3f1VoL6ZtamKfZtZfYIP3qfc/YUymiTlmFUUVzKPWbjPLcAs4MRSL5UcLzOrB7SgCi8d7i0ud9/o7nvCxYeBQ6ognJHAqWaWDUwFjjWzJ0u1ScbxqjCuJB2v4n2vDf9dD7wIDC/VJLa/yVRMCtOAi8O794cBW919XbKDMrMOxddRzWw4we9mYxXs14C/AZ+6+1/30qzKj1mUuJJxzMysrZm1DB83AkYDn5VqNg24JHx8NvCmh3cHkxlXqWvOpxLcp4mVu9/g7l3cPY3gJvKb7n5hqWZVfryixJWM4xXut4mZNSt+DJwAlO61GNvfZL39sZHqxMyeJuiV0sbMcoBbCG664e4PAq8S3LnPAnYBl1aTuM4GfmZmBcBu4Ly4/zBCI4GLgEXh9WiAG4FuCbEl45hFiSsZx6wj8JiZ1SVIQs+4+ytmNgHIdPdpBMnsCTPLIvjGe17MMUWN62ozOxUoCOMaWwVxlakaHK8ocSXreLUHXgy/79QD/u7u/zazqyD+v0mNaBYRkRKpePlIRET2QklBRERKKCmIiEgJJQURESmhpCAiIiWUFKRGM7NZZhb7XLpmdrUFFVufintfyWRmLc3s58mOQ5JHSUFSVjh6NqqfAye7+wVxxVNNtCR4r5KilBQkdmaWFn7LftiCWv8zwlG33/qmb2ZtwrIDxbXsXzKzl81slZmNN7PrLKh9/6GZtUrYxYVm9r6ZLQ5HNhePCp1iZvPCdU5L2O6zZvYyMKOMWK8Lt7PYzK4Nn3sQOAiYZma/KtW+rpndbkHt+4Vm9svw+ePC/S4K42gYPp9tZn8ysw/MLNPMhpnZdDNbUTw4yYJ5It4xsxfNbKmZPWhmdcLXxoTbXGxmf06IY4eZ/dGC4oAfmln78Pm2ZvZ8eBzmmdnI8PnfhXHNMrOVZnZ1uKlbgZ4W1PH/i5l1DGP5JNznD/b5P4LUDPurBrd+9LO3HyCNYFTokHD5GeDC8PEsICN83AbIDh+PJRit2QxoS1A586rwtTsJCuQVr/9w+PhIwvkqgD8l7KMlwXwMTcLt5gCtyojzEGBR2K4psAQYGr6WTRn17YGfEdRnqhcutwIOIKhg2Tt87vGEeLOBnyW8j4UJ73F9+PzRQC5BIqoLvE4wersT8EXYth7wJnB6uI4DPw4f3wbcFD7+OzAqfNyNoGwIwO+A94GG4XHfSDDCPo2EOT+AXwO/DR/XBZol+/+TfuL9qXVlLqTaWuXuxeUq5hN8+FTkLQ/mUthuZlsJqlZC8ME9KKHd0xDMWWFmzS2oAXQCQcGz34RtDiAskUFQVrqsgmujgBfdfSeAmb0A/AD4uJwYRxNMEFMQxrDJzAaH73d52OYx4BcEpZohqFtT/D6aJrzH3DB2gLnuvjKM4+kwtnxglrtvCJ9/iiARvgTkAcWzh80Hjk+Ir5/9Z3qC5hbW1QH+5UHBtz1mtp6gvEJp84ApFhQnfCnhdyi1lJKCVJU9CY8LgUbh4wL+cxnzgHLWKUpYLuLb/3dL12pxgtLCZ7n7ssQXzGwEsHMvMe7LxC5Wxv4r2k7i+yj9Hovf197e097ku3vxOoUJ26kDHO7uu78VYJAkSv9OvvN5ECbaIwkmo3nCzP7i7o+XE4fUcLqnIMmWzX9KEp+9j9s4F8DMRhFUi9wKTAd+aVZSRXVohO28A5xuZo0tqE55BvBuBevMAK4qvmkd3uv4DEgzs15hm4uAtyv5noabWY/wXsK5wGyCSYaOCu+91AXGRNjuDGB88YKZDamg/XaCy1nF7bsTXNZ6mKBwXZXMaS7JozMFSbbbgWfM7CKCa+T7YrOZvQ80By4Ln/sDweWahWFiyAZOKW8j7v6RmT0KzA2fmuzu5V06gmAqx97hfvIJ7m/ca2aXAs+GyWIe8GAl39MHBDd9BxIkqxfdvcjMbgDeIjhreNXd/1nBdq4G7jOzhQR/7+8AV+2tsbtvNLP3zGwxwRzAi4H/F763HcDFlXwfUsOoSqpINWNmRwO/cfdyk5hIHHT5SERESuhMQURESuhMQURESigpiIhICSUFEREpoaQgIiIllBRERKSEkoKIiJT4/x1i0Sbu5283AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the relationship between cumulative variance explained and number of components\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(np.array(range(len(compressor.explained_variance_ratio_))) + 1, \n",
    "         np.cumsum(compressor.explained_variance_ratio_))\n",
    "plt.xlabel('number of components')\n",
    "plt.ylabel('cumulative explained variance')\n",
    "plt.show(); # top 5 components can already explained 90% of the original variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supervised Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LDA (Linear Discriminant Analysis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared with PCA, LDA is a supervised technique attempts to find a feature subset to maximize class linear-separability, that is, the projected points of the observations with the same class label are as close as possible, while the distances between the centers of difference class labels are as large as possible. LDA can only be applied to classification problems. LDA assumes that classes are normally distributed and have the same covariance matrix.  \n",
    "    \n",
    "Math detail can be accessed at the [official website](https://scikit-learn.org/stable/modules/lda_qda.html#lda-qda) of sklearn. Traditionally, LDA will reduce dimension to (K-1) where K is the number of classes. But in sklearn, it allows further dimension by incorporating PCA into LDA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-30T03:37:18.844214Z",
     "start_time": "2020-03-30T03:37:18.787150Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\n",
    "\n",
    "# classification example\n",
    "# use iris dataset\n",
    "from sklearn.datasets import load_iris\n",
    "iris = load_iris()\n",
    "X, y = iris.data, iris.target\n",
    "\n",
    "# random suffle the dataset\n",
    "# use the first 100 observations as train_set\n",
    "# the rest 50 observations as test_set\n",
    "np.random.seed(1234)\n",
    "idx = np.random.permutation(len(X))\n",
    "X = X[idx]\n",
    "y = y[idx]\n",
    "\n",
    "train_set = X[0:100,:]\n",
    "test_set = X[100:,]\n",
    "train_y = y[0:100]\n",
    "test_y = y[100:,]\n",
    "\n",
    "# we need to standardize the data because LDA assumes normal distribution\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "model = StandardScaler()\n",
    "model.fit(train_set) \n",
    "standardized_train = model.transform(train_set)\n",
    "standardized_test = model.transform(test_set)\n",
    "\n",
    "# start compressing\n",
    "compressor = LDA(n_components=2) # set n_components=2\n",
    "# n_components <= min(n_classes - 1, n_features)\n",
    "\n",
    "compressor.fit(standardized_train, train_y) # fit on trainset\n",
    "transformed_trainset = compressor.transform(standardized_train) # transform trainset, (100,2)\n",
    "transformed_testset = compressor.transform(standardized_test) # transform test set\n",
    "assert transformed_trainset.shape[1] == transformed_testset.shape[1] # same number of features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-30T03:37:19.129206Z",
     "start_time": "2020-03-30T03:37:18.847470Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VfX9x/HXh7333nsTBMNwI9WKC0Vq3bOuqrX1VxVw1IG46qitq2pRabVWWaKC4AQVVEAlCTvsAMreBDI+vz/OoUYKyYHk5ma8n4/HfXDvGfe+D4H7yVmfr7k7IiIiR6pMvAOIiEjxpkIiIiL5okIiIiL5okIiIiL5okIiIiL5okIiIiL5okIiIiL5okIiIiL5okIiIiL5Ui7eAQpDvXr1vFWrVvGOISJSrMyZM2eju9fPa7lSUUhatWrF7Nmz4x1DRKRYMbOVUZbToS0REckXFRIREckXFRIREckXFRIREckXFRIREcmXmBYSMxtlZuvNLOUQ883M/mpmqWaWZGa9csy7wsyWhI8rckw/2sySw3X+amYWy20QEZHcxXqP5FVgYC7zTwfah4/rgOcBzKwOcC/QF+gD3GtmtcN1ng+X3b9ebu8vIiIxFtNC4u7Tgc25LHIOMNoDXwG1zKwxcBrwobtvdvctwIfAwHBeDXef6cEYwaOBc2O5DSIixdGarXu4/915ZGZlx/yz4n1DYlNgdY7XaeG03KanHWT6/zCz6wj2XGjRokXBJRYRKcKys53Xv17JI5MXku0wuGdTEprViulnxruQHOz8hh/B9P+d6P4i8CJAYmLiQZcRESlJlm3YybCxyXyzYjMntK/HQ4O707xOlZh/brwLSRrQPMfrZsDacHr/A6Z/Fk5vdpDlRURKrcysbF76fDlPfbSYSuXK8OdfJfCro5tRWNcixfvy34nA5eHVW/2Abe6+DpgC/NLMaocn2X8JTAnn7TCzfuHVWpcD78QtvYhInM1bu41zn/uSRz9YyICODfjojydxfmLzQisiEOM9EjP7N8GeRT0zSyO4Eqs8gLu/AEwCzgBSgd3AVeG8zWY2ApgVvtUD7r7/pP1vCa4GqwxMDh8iIqVKekYWf/tkCS9MW0btKhV4/pJenN69cVyyWHDxU8mWmJjo6v4rIiXFnJWbuWNMEks37GJIr2bcc1ZnalWpUOCfY2Zz3D0xr+XifY5EREQi2rU3kz9PWcRrM1fQpGZlXru6Dyd1yHO4kJhTIRERKQamL97A8HHJrN22hyuOacVtp3WkWsWi8RVeNFKIiMhBbd29jwffX8CYOWm0qV+Vt68/hsRWdeId62dUSEREiqjJyeu45515bNm9j5tObsvvBrSnUvmy8Y71P1RIRESKmPU70rn3nXlMTvmBrk1q8NrVvenapGa8Yx2SComISBHh7oyZk8aD7y9gT0YWdwzsyLUntKF82Xjf8pc7FRIRkSJg9ebd3Dk+mc+XbKR3q9o8MiSBtvWrxTtWJCokIiJxlJ3tjJ65gsemLMKAEed05ZK+LSlTpvgMtaRCIiISJ6nrdzB0bDJzVm7hpA71GTm4G81qx77JYkFTIRERKWQZWdm8OH0ZT3+0hCoVy/Lkr3swuGfTQu2PVZBUSEREClHKmm3cMSaJ+eu2c2b3xtw3qCv1q1eMd6x8USERESkE6RlZPP3xEl6cvow6VSvwwqVHM7Bbo3jHKhAqJCIiMTZrxWaGjkli2cZdXJDYnDvP6EzNKuXjHavAqJCIiMTIzr2ZPPbBQkbPXEmz2pX512/6cnz7evGOVeBUSEREYuDTReu5a1wy67anc/VxrbnttA5UqVAyv3JL5laJiMTJll37GPHefMZ9t4Z2Daox5oZjObpl7XjHiikVEhGRAuDuTEr+gXsnprB1dwa3DGjHTQPaUbFc0WuyWNAiFRIzOx5o7+6vmFl9oJq7L49tNBGR4mH99nTunpDC1Pk/0r1pTUZf3ZcuTWrEO1ahybOQmNm9QCLQEXiFYMz1fwHHxTaaiEjR5u68PTuNEe/PZ19mNsNP78Rvjm9NuSLeZLGgRdkjGQz0BL4FcPe1ZlY9pqlERIq4VZt2M3x8El+mbqJP6zo8OiSB1vWqxjtWXEQpJPvc3c3MAcysdP5NiYgAWdnOqzNW8PiURZQtYzx4bjcu7tOiWDVZLGhRCslbZvZ3oJaZXQtcDbwU21giIkXPkh93cMfYJL5btZWTO9Zn5ODuNKlVOd6x4i7PQuLuj5vZqcB2gvMkf3L3D2OeTESkiNiXmc0L05byzCepVK1Ylr9ccBTnHNWk2DZZLGhRTra3Bj7fXzzMrLKZtXL3FbEOJyISb0lpW7ljTBILf9jB2T2acO/ZXahXrXg3WSxoUQ5tvQ0cm+N1Vjitd0wSiYgUAXv2ZfGXjxbz0ufLqF+9Ii9dnsipXRrGO1aRFKWQlHP3fftfuPs+M6sQw0wiInH11bJNDBubxIpNu7moT3OGn9GZGpVKTpPFghalkGwws0HuPhHAzM4BNsY2lohI4duRnsEjkxfy+teraFGnCm9c05dj25W8JosFLUohuQF43cyeAQxYDVwe01QiIoXsk4U/ctf4FH7cns41x7fmj7/sSOUKJb+9SUGIctXWUqCfmVUDzN13xD6WiEjh2LxrHw+8O48J36+lQ8NqPHfJsfRsUbKbLBa0KFdtVQSGAK2Acvsvd3P3B2KaTEQkhtydd5PWcd/EeexIz+APp7Tnxv7tqFCudLU3KQhRDm29A2wD5gB7YxtHRCT2ftiWzt0TkvlowXp6NK/FY0MS6NhInZ+OVJRC0szdB8Y8iYhIjLk7b85azUPvLyAjO5u7z+zMVce1pmwpbm9SEKLsw80ws+5H8uZmNtDMFplZqpkNO8j8lmb2sZklmdlnZtYsx7xHzSwlfFyQY/ovzOxbM/vezL4ws3ZHkk1ESpeVm3Zx8UtfM3xcMt2a1mTKH07kmhPaqIgUgCh7JMcDV5rZcoJDWwa4uyfktpKZlQWeBU4F0oBZZjbR3efnWOxxYLS7v2ZmA4CHgcvM7EygF3AUUBGYZmaT3X078DxwjrsvMLMbgbuBK6NvsoiUJlnZzitfLufxqYsoX6YMD5/XnQt7N1d7kwIUpZCcfoTv3QdIdfdlAGb2JnAOkLOQdAFuDZ9/CkzIMX2au2cCmWY2FxgIvAU4sH/EmJrA2iPMJyIl3KIfgiaLc1dv5ZTODXjw3O40qlkp3rFKnCiX/64EMLMGwOH8BJoS3HOyXxrQ94Bl5hJcEfY0wbgn1c2sbjj9XjN7EqgCnMxPBegaYJKZ7SFoJNnvMDKJSCmwLzObZz9N5bnPUqleqTx/vagnZyc01l5IjES5/HcQ8ATQBFgPtAQWAF3zWvUg0/yA17cBz5jZlcB0YA2Q6e5Tzaw3MAPYAMwEMsN1bgXOcPevzex24EmC4nJg7uuA6wBatGiRR1QRKSm+X72VO8bMZfGPOzn3qCb86eyu1Kmqrk6xFOXQ1giC3/o/cveeZnYycFGE9dKA5jleN+OAw1DuvhY4DyC84XGIu28L540ERobz3gCWhOPF93D3r8O3+A/wwcE+3N1fBF4ESExMPLCAiUgJs2dfFk9MXcSoL5fTsEYlRl2ZyIBOarJYGKIUkgx332RmZcysjLt/amaPRlhvFtA+bEO/BrgQuDjnAmZWD9js7tnAcGBUOL0sUCv83AQgAZgarlbTzDq4+2KCE/kLImQRkRJsxtKNDBubzKrNu7mkbwuGnd6J6mqyWGiiFJKt4d7CdIKeW+v56TDTIbl7ppndDEwBygKj3H2emT0AzA6bQPYHHg6H8Z0O3BSuXh74PDyeuR24NDzxTjhK41gzywa2EIzYKCKl0Pb0DB6etIB/f7OaVnWr8OZ1/ejXpm68Y5U65p77UZ9wjPZ0gnMelxBcKfW6u2+KfbyCkZiY6LNnz453DBEpQB/O/5G7JySzYcderj2hDX84pYOaLBYwM5vj7ol5LRflqq1dOV6+lq9UIiL5tHHnXu6bOI/3ktbRqVF1Xro8kYRmteIdq1Q7ZCExsy/c/Xgz28HPr7baf0NijUOsKiJS4Nydd75fy/3vzmPX3iz+eGoHrj+prZosFgGHLCTufnz4pzqZiUhcrd26h7snpPDJwvX0bBE0WWzfUF9NRUWuh7bMrAyQ5O7dCimPiMh/ZWc7b3yzikcmLyQr2/nTWV244thW6o9VxORaSNw928zmmlkLd19VWKFERJZv3MWwsUl8vXwzx7Wry8ODE2hRt0q8Y8lBRLn8tzEwz8y+Af574t3dB8UslYiUWplZ2fzji+U8+eFiKpQrw2NDEjg/sZnamxRhUQrJ/TFPISICzF+7naFjk0hes41fdmnIiHO70bCGmiwWdVEu/51WGEFEpPTam5nFM5+k8vxnS6lVpTzPXtyLM7o30l5IMRGlaWM/4G9AZ6ACwV3qu3T5r4gUhDkrtzB0bBKp63dyXq+m3HNmF2qryWKxEuXQ1jMEfbLeBhKBy4H2sQwlIiXf7n2Z/HnKIl6dsYLGNSrxylW9Obljg3jHkiMQpZDg7qlmVtbds4BXzGxGjHOJSAn2xZKNDBuXRNqWPVx+TEvuGNiJahUjfR1JERTlJ7fbzCoA35vZY8A6oGpsY4lISbRtdwYjJ83nrdlptKlXlbeuP4Y+revEO5bkU5RCchlQBriZYFCp5gSjGoqIRPZByg/c804Km3ft47f92/L7X7SnUnk1WSwJohSSXsAkd9+OLgUWkcO0YUfQZPH95HV0aVyDV67sTbemNeMdSwpQlEIyCPiLmU0H3gSm7B8bRETkUNydcd+u4YH35rNnXxa3n9aR605sQ/myarJY0kS5j+QqMysPnE4wwuFzZvahu//POOkiIgBrtu7hznHJTFu8gaNb1ubRIQm0a1At3rEkRqJetZVhZpMJ2slXBs4BVEhE5Geys51/fb2SRycvxIH7zu7C5ce0ooyaLJZoUW5IHEhwH8nJwGfAy8CvYxtLRIqbpRt2MmxsErNWbOGE9vV4aHB3mtdRk8XSIMoeyZUE50aud/e9sY0jIsVNRlY2L32+jL98tITK5cvy+Pk9GNKrqdqblCJRzpFcWBhBRKT4SVmzjaFjk5i3djund2vE/ed0pUF1NVksbXQrqYgctvSMLP72yRJemLaM2lUq8PwlvTi9e+N4x5I4USERkcMye8Vm7hibxLINu/jV0c24+8zO1KqiJoulmQqJiESya2/QZPG1mStoUrMyo6/uw4kd6sc7lhQBhywkZpZMcLnvQbl7QkwSiUiRM23xBu4cl8zabXu44phW3H5aR6qqyaKEcvuXcFb4503hn/8M/7wE2B2zRCJSZGzdvY8R7y1g7LdptK1flbevP4bEVmqyKD93yELi7isBzOw4dz8ux6xhZvYl8ECsw4lI/ExOXsc978xjy+593HxyO24e0E5NFuWgouybVjWz4939CwAzOxa1kRcpsdZvT+dP78zjg3k/0LVJDV67ujddm6jJohxalELyG2CUmdUkOGeyDbg6pqlEpNC5O2PmpDHivfmkZ2YzdGAnrj2hNeXUZFHyEOWGxDlADzOrAZi7b4t9LBEpTKs37+bO8cl8vmQjfVrV4eEh3WlbX00WJZoovbYaAg8BTdz9dDPrAhzj7v+IeToRiamsbGf0zBX8ecoiDBhxTlcu6dtSTRblsEQ5tPUq8ApwV/h6MfAfQIVEpBhLXb+DoWOTmbNyCyd1qM9D53Wnaa3K8Y4lxVCUQlLP3d8ys+EA7p5pZlkxziUiMZKRlc3fpy3lrx+nUqViWZ78dQ8G91STRTlyUQrJLjOrS3hzopn1IzjhLiLFTMqabdw+JokF67ZzZkJj7ju7K/WrV4x3LCnmolyO8X/ARKBteP/IaOB3Ud7czAaa2SIzSzWzYQeZ39LMPjazJDP7zMya5Zj3qJmlhI8Lckw3MxtpZovNbIGZ3RIli0hplp6RxSOTF3LOs1+ycede/n7Z0Tx7cS8VESkQUa7a+tbMTgI6AgYscveMvNYzs7LAs8CpQBowy8wmuvv8HIs9Dox299fMbADwMHCZmZ0J9AKOAioC08xssrtvJxgfpTnQyd2zzazBYWyvSKnz9bJNDBuXzPKNu7ggsTl3ntGZmlXKxzuWlCBRm+X0AVqFy/cyM9x9dIR1Ut19GYCZvUkwRG/OQtIFuDV8/ikwIcf0ae6eCWSa2VxgIPAW8FvgYnfPBnD39RG3QaRU2ZGewWMfLOKfX62keZ3KvH5NX45rVy/esaQEinL57z+BtsD3wP6T7E5wiCs3TYHVOV6nAX0PWGYuMAR4GhgMVA/Px8wF7jWzJ4EqBMP87i9AbYELzGwwsAG4xd2X5LUdIqXJp4vWc9e4ZNZtT+fq41pz22kdqFJBTRYlNqL8y0oEurj7ITsBH8LBLgE58D1uA54xsyuB6cAaINPdp5pZb2AGQbGYCWSG61QE0t090czOA0YBJ/zPh5tdB1wH0KJFi8OMLlI8bdm1jxHvzWfcd2to36AaY397LL1a1I53LCnhohSSFKARsO4w3zuN4FzGfs2AtTkXcPe1wHkAZlYNGLL/znl3HwmMDOe9Aezf60gDxobPxxPc4/I/3P1F4EWAxMTEwy2CIsWKu/N+8jrufWce2/ZkcMuAdtw0oB0Vy6nJosRepPtIgPlm9g2wd/9Edx+Ux3qzgPZm1ppgT+NC4OKcC5hZPWBzeL5jOMHexf4T9bXcfZOZJQAJwNRwtQnAgHDZkwhukBQptX7cns7dE1L4cP6PdG9ak39d05fOjWvEO5aUIlEKyX1H8sbhjYs3A1OAssAod59nZg8As919ItAfeNjMnODQ1v6xT8oDn4c3SG0HLg1PvAM8ArxuZrcCO4FrjiSfSHHn7rw1ezUPvr+AfZnZ3HlGJ64+Tk0WpfDZ4Z/6KH4SExN99uzZ8Y4hUmBWbdrNsHFJzFi6ib6t6/DokARa1dPoDlKwzGyOuyfmtVxuQ+1+4e7Hm9kOfn6S3AB3d+07ixSyrGzn1RkreHzKIsqWMUYO7sZFvVuoyaLEVW4jJB4f/lm98OKIyKEs/nEHd4xJ4vvVWxnQqQEjB3ejcU01WZT4i3xheXgHeaX9r919VUwSicjP7MvM5vnPlvLMp0uoVrEcT194FIN6NFGTRSkyotyQOAh4AmgCrAdaAguArrGNJiJzV29l6NgkFv6wg7N7NOG+s7tQt5r6Y0nREmWPZATQD/jI3Xua2cnARbGNJVK67dmXxVMfLeblz5dRv3pFXro8kVO7NIx3LJGDilJIMsL7OcqYWRl3/9TMHo15MpFSaubSTQwfl8SKTbu5qE8Lhp/RiRqV1GRRiq4ohWRreNf5dIL7N9bzU7sSESkg29MzeGTyQt74ehUt61bhjWv7cmxbNVmUoi9KITkHSCfo0nsJUBN4IJahREqbTxb+yJ3jUli/I51rT2jN/53akcoV1N5Eioco45HsyvHytRhmESl1Nu3cywPvzeed79fSsWF1XrjsaI5qXivesUQOS243JB70RkR0Q6JIvrk7E+eu5f5357MjPYM/nNKeG/u3o0I5tTeR4ie3GxJ1I6JIDKzbtoe7x6fw8cL19Ghei8eGJNCxkf67SfEV6YZEM+sFHE+wR/KFu38X01QiJVB2tvPmrNU8PGkBGdnZ3H1mZ646rjVl1d5EirkoNyT+CTgfGBdOetXM3nb3B2OaTKQEWbFxF8PGJfHVss0c06YujwzpTsu6arIoJUOUPZKLgJ7ung5gZo8A3wIqJCJ5yMzK5pUvV/DEh4soX6YMj5zXnQt6N1d7EylRohSSFQQ9ttLD1xWBpbEKJFJSLPxhO0PHJDE3bRundG7Ag+d2p1HNSnmvKFLMRCkke4F5ZvYhwTmSU4EvzOyvAO5+SwzziRQ7ezOzePbTpTz3aSo1K5fnbxf15KyExtoLkRIrSiEZHz72+yw2UUSKv+9WbWHo2CQW/7iTwT2bcs9ZXahTtUK8Y4nEVJRCMtnd1+ecYGYd3X1RjDKJFDu792XyxNTFjPpyOY1qVGLUlYkM6KQmi1I6RCkkn5vZPe7+FoCZ/RH4DdAlpslEiokZqRsZNi6ZVZt3c2m/Fgwd2InqarIopUiUQtIfeNHMzgcaEoxF0ieWoUSKg217Mnh40gLenLWaVnWr8OZ1/ejXpm68Y4kUuii9ttaZ2QfAcCAbGO7uO2OeTKQImzrvB+6ekMLGnXu5/qQ23HpKByqVV5NFKZ2i3JD4IbAO6AY0A0aZ2XR3vy3W4USKmo0793LfxHm8l7SOTo2q8/IViSQ0U5NFKd2iHNp61t0nhM+3mtmxBHsnIqWGuzPh+zXc/+58du/N4o+nduCG/m0pX1ZNFkWiHNqaYGYtgfbu/hFQHvhLzJOJFBFrt+7hrvHJfLpoAz1bBE0W2zdUk0WR/aIc2roWuA6oA7QlOLz1AvCL2EYTia/sbOf1b1bx6OSFZGU7fzqrC1cc20pNFkUOEOXQ1k0EV2l9DeDuS8ysQUxTicTZsg07GTY2mW9WbOb4dvV4+LzuNK9TJd6xRIqkSC1S3H3f/vYOZlaOnw94JVJiZGZl8/IXy3nqw8VULFeGx36VwPlHN1N7E5FcRCkk08zsTqCymZ0K3Ai8G9tYIoVv/trt3DF2LilrtnNa14aMOKcbDWqoyaJIXqIUkmEEd7InA9cDk4CXYxlKpDDtzczimU9Sef6zpdSqUp7nLunF6d0aaS9EJKIoV21lAy+FD5ESZc7KoMli6vqdnNerKfec2YXaarIoclgiDbUrUtLs2pvJ41MX8eqMFTSpWZlXr+pN/466hkTkSKiQSKnz+ZINDB+XTNqWPVx+TEvuGNiJahX1X0HkSEX+32NmVd19VyzDiMTStt0ZPPj+fN6ek0abelV56/pj6NO6TrxjiRR7efZ3MLNjzWw+QddfzKyHmT0X5c3NbKCZLTKzVDMbdpD5Lc3sYzNLMrPPzKxZjnmPmllK+LjgIOv+zczUPFIi+SDlB055ahrjvlvDjf3bMun3J6iIiBSQKHskTwGnARMB3H2umZ2Y10pmVhZ4lmBo3jRglplNdPf5ORZ7HBjt7q+Z2QDgYeAyMzsT6AUcRTBG/DQzm+zu28P3TgTUKU/ytH5HOvdNnMek5B/o0rgGr1zZm25Na8Y7lkiJEunQlruvPuBSyKwIq/UBUt19GYCZvQmcA+QsJF2AW8PnnwITckyf5u6ZQKaZzQUGAm+FBerPwMXA4Cj5pfRxd8Z+u4YR781nT0YWt5/WketObKMmiyIxEOV/1eqw46+bWQUzu43wMFcemgKrc7xOC6flNBcYEj4fDFQ3s7rh9NPNrIqZ1QNOBpqHy90MTHT3dREySCmUtmU3V7wyi9venku7BtWYdMsJ3HRyOxURkRiJskdyA/A0QRFIA6YS9N/Ky8Hu5jqwtcptwDNmdiUwHVgDZLr7VDPrDcwANgAzCfZMmgDnE4zamPuHm11H0GySFi1aRIgrxV12tvPPr1by6AcLAbh/UFcu69eSMmqyKBJTUQqJufslR/Deafy0FwFB1+C1ORdw97XAeQBmVg0Y4u7bwnkjgZHhvDeAJUBPoB2QGh5qq2Jmqe7e7sAPd/cXgRcBEhMT1RushFu6YSdDxyQxe+UWTuxQn4cGd6NZbTVZFCkMUQrJDDNbDvwHGOvuWyO+9yygvZm1JtjTuJDgvMZ/hYetNod3zw8HRoXTywK13H2TmSUACcDU8JxJoxzr7zxYEZHSIyMrmxenL+Ppj5dQuXxZHj+/B0N6NVV7E5FCFKVFSnsz60NQCO4KLwV+093/lcd6mWZ2MzAFKAuMcvd5ZvYAMNvdJxIconrYzJzg0Nb+Q2blgc/DL4PtwKVhERH5r5Q12xg6Nol5a7dzRvdG3DeoKw2qq8miSGEz9+hHfcI9iCeBS9y9bMxSFbDExESfPXt2vGNIAUnPyOKvHy/h79OXUbtKBR48tysDuzWOdyyREsfM5rh7Yl7LRRkhsQbBFVUXEoyQOJ7g0l6RQjdrxWaGjkli2cZdnH90M+4+sws1q5SPdyyRUi3KOZK5BPd3PODuM2OcR+Sgdu7N5LEPFjJ65kqa1qrM6Kv7cGKH+vGOJSJEKyRt/HCOf4kUsGmLN3DnuGTWbtvDlce24vbTOlJVTRZFioxD/m80s7+4+x+AieHJ8J9x90ExTSal3tbd+3jgvfmM+3YNbetXZcwNx3B0S/XHEilqcvu17p/hn48XRhCRnCYlr+NP76SwdXcGN5/cjpsHtKNS+WJzfYdIqXLIQuLuc8KnR7n70znnmdnvgWmxDCal0/rt6dzzTgpT5v1It6Y1eO3qPnRtoiaLIkVZlAPNVxC0SMnpyoNMEzli7s7bc9J48L35pGdmM3RgJ649oTXl1B9LpMjL7RzJRQR3orc2s4k5ZlUHNsU6mJQeqzfvZvi4ZL5I3UifVnV4ZEh32tSvFu9YIhJRbnskM4B1QD3giRzTdwBJsQwlpUNWtjN65goe+2ARZQxGnNuNS/q0UJNFkWImt3MkK4GVwDGFF0dKi9T1O7hjTBLfrtpK/471GTm4O01rVY53LBE5AlHubO8H/A3oDFQg6Ju1y91rxDiblEAZWdm88NlS/vZJKlUqluWpC3pw7lFqsihSnEU52f4MQXuUt4FE4HKCVu4ihyU5bRu3j5nLwh92cGZCY+4f1JV61SrGO5aI5FPUoXZTzaysu2cBr5jZjBjnkhIkPSOLpz5azEvTl1GvWkX+ftnRnNa1Ud4rikixEKWQ7DazCsD3ZvYYwQn4qrGNJSXF18s2MWxcMss37uLC3s0ZfkZnalZWk0WRkiRKIbmM4LzIzcCtBKMeDsl1DSn1dqRn8OgHC/nXV6toXqcyr1/Tl+Pa1Yt3LBGJgSgDW60Mn+4B7o9tHCkJPl24nrvGJ7Nuezq/Ob41f/xlB6pUUJNFkZIqtxsSk4FDdv1194SYJJJia/OufYx4bz7jv1tD+wbVGPvbY+nVona8Y4lIjOX2a+JZhZZCijV3572kddw3cR7b9mRwyy/ac9PJbalYTk0WRUqDvG5IFMnVj9vTuWt8Ch8t+JGEZjX51zV96dxYtxiJlCZRbkjcwU+HuCoA5dENiaWeu/MJuwFoAAAQyUlEQVSfWasZOWkB+zKzufOMTlx9nJosipRGUU62V8/52szORWO2l2qrNu1m2LgkZizdRN/WdXh0SAKt6umKcJHS6rAvpXH3CWY2LBZhpGjLynZe+XI5j09dRLkyZXhocHcu7N1cTRZFSrkoh7bOy/GyDEGbFI3hXsos+mEHd4xNYu7qrQzo1ICRg7vRuKaaLIpItD2Ss3M8zwRWAOfEJI0UOfsys3nus1Se/TSV6pXK8/SFRzGoRxM1WRSR/4pyjuSqwggiRc/c1Vu5Y0wSi37cwaAeTbj37C7UVZNFETlAlENbrYHfAa1yLu/ug2IXS+Jpz74snvxwEf/4YjkNqlfi5csTOaVLw3jHEpEiKsqhrQnAP4B3gezYxpF4m7l0E8PGJbFy024u7tuCYad3okYlNVkUkUOLUkjS3f2vMU8icbU9PYOHJy3k39+somXdKrxxbV+ObasmiyKStyiF5GkzuxeYCuzdP9Hdv41ZKilUH83/kbsmJLNhx16uO7ENt57SgcoV1N5ERKKJUki6E7SSH8BPh7Y8fC3F2Kade7n/3flMnLuWjg2r8/fLEjmqea14xxKRYiZKIRkMtHH3fbEOI4XD3Zk4dy33TZzHzr2Z3HpKB37bvy0Vyqm9iYgcviiFZC5QC1gf4yxSCNZt28Pd41P4eOF6jmpei8d+lUCHhtXzXlFE5BCiFJKGwEIzm8XPz5Ho8t9iJDvb+fesVTw8aSGZ2dncfWZnrjquNWXV3kRE8ilKIbn3SN/czAYCTxMM1fuyuz9ywPyWwCigPrAZuNTd08J5jwJnhouOcPf/hNNfJ2jTkgF8A1zv7hlHmrE0WL5xF8PGJvH18s0c27Yuj5yXQIu6VeIdS0RKiCh3tk87kjc2s7LAs8CpQBowy8wmuvv8HIs9Dox299fMbADwMHCZmZ0J9AKOAioC08xssrtvB14HLg3XfwO4Bnj+SDKWdJlZ2Yz6cjlPTF1MhbJleOS87lzQu7nam4hIgYrleCR9gFR3Xxa+z5sEPbpyFpIuwK3h808Jbn7cP32au2cCmWY2FxgIvOXuk3Jk+wZoltc2lEYL1m1n6NgkktK2cUrnhjx4bjca1awU71giUgLleZmOu1d39xrhoxIwBHgmwns3BVbneJ0WTstpbvh+EFwdVt3M6obTTzezKmZWDzgZaJ5zRTMrT3BZ8gcRspQaezOzePLDxZz9ty9Ys2UPz1zck5cuP1pFRERiJpbjkRzs+MmB7edvA54xsyuB6cAaINPdp5pZb2AGsAGYSdB5OKfngOnu/vlBP9zsOuA6gBYtWkSIW/x9u2oLQ8cksWT9Tgb3bMqfzupC7aoV4h1LREq4WI5HksbP9yKaAWtzLuDua4Hzws+pBgxx923hvJHAyHDeG8CSHJnuJThBf/2hPtzdXwReBEhMTCzR46fs3pfJE1MXM+rL5TSqUYlXruzNyZ0axDuWiJQSsRyPZBbQPuwevAa4ELg45wLhYavN7p4NDCe4gmv/ifpa7r7JzBKABIIWLZjZNcBpwC/C9Uq1L1M3MmxcEqs37+HSfi0YOrAT1dVkUUQKUczGI3H3TDO7GZhCcPnvKHefZ2YPALPdfSLQH3jYzJzg0NZN4erlgc/Dq4u2E1wWvP/Q1gvASmBmOH+cuz9wJBmLs217Mnjo/QX8Z/ZqWteryn+u60ffNnXjHUtESiFzz/2oj5m9Bvze3beGr2sDT7j71YWQr0AkJib67Nmz4x2jwEyd9wN3T0hh4869XBs2WaxUXk0WRaRgmdkcd0/Ma7koh7YS9hcRAHffYmY985VOjsiGHXu57915vJ+0jk6NqvPyFYkkNFOTRRGJryiFpIyZ1Xb3LQBmVifielJA3J0J36/h/nfns3tvFrf9sgPXn9SW8mXVZFFE4i9KQXgCmGFmYwiu1vo14dVUEntrtu7hrvHJfLZoA71aBE0W2zVQk0URKTqinGwfbWazCcYfMeC8A9qcSAxkZzuvf72SRyYvJNvh3rO7cPkxrdRkUUSKnEiHqMLCoeJRSJZt2Mmwscl8s2Izx7erx8Pndad5HTVZFJGiSec6ipDMrGxe+nw5T320mErlyvDYrxI4/+hmarIoIkWaCkkRMX/tdu4YO5eUNds5rWtDRpzTjQY11B9LRIo+FZI4S8/I4plPUnlh2lJqVanA85f04vTujeMdS0QkMhWSOJqzcjN3jEli6YZdDOnVjHvO6kytKmqyKCLFiwpJHOzam8mfpyzitZkraFKzMq9d3YeTOtSPdywRkSOiQlLIpi/ewPBxyazZuocrjmnJ7QM7Ua2ifgwiUnzpG6yQbNudwYj35zNmThpt6lfl7RuOoXerOvGOJSKSbyokheCDlHXc8848Nu/ax43923LLL9qryaKIlBgqJDG0fkc6974zj8kpP9ClcQ1eubI33ZrWjHcsEZECpUISA+7OmDlpPPj+AvZkZHH7aR257sQ2arIoIiWSCkkBW715N3eOT+bzJRtJbFmbR4Yk0K5BtXjHEhGJGRWSApKd7YyeuYLHpiwC4P5BXbmsX0vKqMmiiJRwKiQFIHX9ToaNTWL2yi2c2KE+Dw3uRrPaarIoIqWDCkk+ZGRl8+L0ZTz90RIqVyjLE+f34LxeTdVkUURKFRWSI5SyZht3jEli/rrtnNG9EfcP6kb96hXjHUtEpNCpkBym9Iwsnv54CS9OX0adqhV44dJeDOymJosiUnqpkByGWSs2M3RMEss27uL8o5tx95ldqFmlfLxjiYjElQpJBDv3ZvLYBwsZPXMlzWpX5p+/6cMJ7dVkUUQEVEjy9Nmi9dw1PoW12/Zw1XGtuO2XHamqJosiIv+lb8RcDB+XzL+/WUW7BtUYc8OxHN2ydrwjiYgUOSokuWhVtwq/G9COmwe0o2I5NVkUETkYFZJcXH9S23hHEBEp8tRFUERE8kWFRERE8kWFRERE8kWFRERE8kWFRERE8kWFRERE8kWFRERE8kWFRERE8sXcPd4ZYs7MNgArj3D1esDGAoxTHGibSwdtc8mX3+1t6e55dqgtFYUkP8xstrsnxjtHYdI2lw7a5pKvsLZXh7ZERCRfVEhERCRfVEjy9mK8A8SBtrl00DaXfIWyvTpHIiIi+aI9EhERyRcVEsDMRpnZejNLOcR8M7O/mlmqmSWZWa/CzljQImzzJeG2JpnZDDPrUdgZC1pe25xjud5mlmVmvyqsbLESZZvNrL+ZfW9m88xsWmHmi4UI/7Zrmtm7ZjY33OarCjtjQTKz5mb2qZktCLfn9wdZJqbfYSokgVeBgbnMPx1oHz6uA54vhEyx9iq5b/Ny4CR3TwBGUDKOLb9K7tuMmZUFHgWmFEagQvAquWyzmdUCngMGuXtX4PxCyhVLr5L7z/kmYL679wD6A0+YWYVCyBUrmcAf3b0z0A+4ycy6HLBMTL/DVEgAd58ObM5lkXOA0R74CqhlZo0LJ11s5LXN7j7D3beEL78CmhVKsBiK8HMG+B0wFlgf+0SxF2GbLwbGufuqcPliv90RttmB6mZmQLVw2czCyBYL7r7O3b8Nn+8AFgBND1gspt9hKiTRNAVW53idxv/+oEqy3wCT4x0i1sysKTAYeCHeWQpRB6C2mX1mZnPM7PJ4ByoEzwCdgbVAMvB7d8+Ob6SCYWatgJ7A1wfMiul3mMZsj8YOMq1UXO5mZicTFJLj452lEPwFGOruWcEvq6VCOeBo4BdAZWCmmX3l7ovjGyumTgO+BwYAbYEPzexzd98e31j5Y2bVCPam/3CQbYnpd5gKSTRpQPMcr5sR/DZToplZAvAycLq7b4p3nkKQCLwZFpF6wBlmlunuE+IbK6bSgI3uvgvYZWbTgR5ASS4kVwGPeHDvQ6qZLQc6Ad/EN9aRM7PyBEXkdXcfd5BFYvodpkNb0UwELg+vfOgHbHP3dfEOFUtm1gIYB1xWwn87/S93b+3urdy9FTAGuLGEFxGAd4ATzKycmVUB+hIcYy/JVhHsgWFmDYGOwLK4JsqH8FzPP4AF7v7kIRaL6XeY9kgAM/s3wdUb9cwsDbgXKA/g7i8Ak4AzgFRgN8FvNMVahG3+E1AXeC78DT2zuDe7i7DNJU5e2+zuC8zsAyAJyAZedvdcL48u6iL8nEcAr5pZMsEhn6HuXpw7Ah8HXAYkm9n34bQ7gRZQON9hurNdRETyRYe2REQkX1RIREQkX1RIREQkX1RIREQkX1RIREQkX1RIpNQJ24HE/FJmM7sl7Mj6eqw/K57MrJaZ3RjvHBI/KiQih8HMDufeqxuBM9z9kljlKSJqEWyrlFIqJFIkmVmr8Lf5l8IxFqaaWeVw3n/3KMysnpmtCJ9faWYTwrEmlpvZzWb2f2b2nZl9ZWZ1cnzEpeE4Kylm1idcv2o4lsWscJ1zcrzv22b2LjD1IFn/L3yfFDP7QzjtBaANMNHMbj1g+bJm9riZJYdjQ/wunP6L8HOTwxwVw+krzOwhM5tpZrPNrJeZTTGzpWZ2Q7hMfzObbmbjzWy+mb1gZmXCeReF75liZo/myLHTzEZaMC7HV+Fd3phZfTMbG/49zDKz48Lp94W5PjOzZWZ2S/hWjwBtLRjT5M9m1jjM8n34mScc8T8EKR7cXQ89itwDaEXQ2vuo8PVbwKXh88+AxPB5PWBF+PxKgjt3qwP1gW3ADeG8pwia2e1f/6Xw+YlASvj8oRyfUYug31TV8H3TgDoHyXk0QQfZqgQtyecBPcN5K4B6B1nntwR9kcqFr+sAlQi6s3YIp43OkXcF8Nsc25GUYxvXh9P7A+kExass8CHwK6AJQUuQ+gSdLD4Bzg3XceDs8PljwN3h8zeA48PnLQhabwDcB8wAKoZ/75sI7hhvtf/vMFzuj8Bd4fOyQPV4/3vSI7YPtUiRomy5u+9v+TCH4AsrL596MCbDDjPbBrwbTk8GEnIs928Ixq4wsxoWDPD0S2CQmd0WLlOJsM0E8KG7H2yMi+OB8R40PcTMxgEnAN/lkvEU4AV3zwwzbLZgBMrl/lNfs9cIBmD6S/h6Yo7tqJZjG9PD7ADfuPuyMMe/w2wZwGfuviGc/jpB8ZwA7APeC9edA5yaI18X+6kDcg0zqx4+f9/d9wJ7zWw90PAg2zcLGGVBI8EJOX6GUkKpkEhRtjfH8yyCNucQ7KnsPyxbKZd1snO8zubn/94P7A3kBH2Xhrj7opwzzKwvsOsQGY+k37wd5PPzep+c23HgNu7frkNt06FkuPv+dbJyvE8Z4Bh33/OzgEFhOfBn8j/fIWFxPhE4E/inmf3Z3UfnkkOKOZ0jkeJoBcEhJQgO3xyJCwDM7HiCTqjbCIbX/Z2F35hm1jPC+0wHzjWzKmZWlWBgrM/zWGcqcMP+E/fhuZuFQCszaxcucxlwuOOn9zGz1uG5kQuALwgGODopPJdUFrgowvtOBW7e/8LMjspj+R0Eh9r2L9+S4JDbSwRdaQt0fHAperRHIsXR48BbZnYZwTH/I7HFzGYANYCrw2kjCA4lJYXFZAVwVm5v4u7fmtmr/DSWxcvuntthLQjGeOkQfk4GwfmaZ8zsKuDtsMDM4vBHapxJcOK7O0GBG+/u2WY2HPiUYO9kkru/k8f73AI8a2ZJBN8R04EbDrWwu28ysy/NLIVgJM0U4PZw23YCpWHUxVJN3X9FSgAz6w/c5u65Fj6RWNChLRERyRftkYiISL5oj0RERPJFhURERPJFhURERPJFhURERPJFhURERPJFhURERPLl/wFuz7e2mCP1vAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the relationship between cumulative variance explained and number of components\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(np.array(range(len(compressor.explained_variance_ratio_))) + 1, \n",
    "         np.cumsum(compressor.explained_variance_ratio_))\n",
    "plt.xlabel('number of components')\n",
    "plt.ylabel('cumulative explained variance')\n",
    "plt.show();  # LDA compresses the original 4 variables into 2 \n",
    "# These 2 variables can explain 100% of the variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "690.488px",
    "left": "28.9922px",
    "top": "134px",
    "width": "319.961px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
